
% Classical Minimum Distance (CMD) gradient

function g = cmd_gradient4(theta,auxdata)

    % Extract vector of sample moments
    mhat        = auxdata{1};
    
    % Extract vector of sample moments from SCF
    mhat_scf    = auxdata{2};
    
    % Extract the fixed income asset share of group 1
    sa1         = auxdata{3};
    
    %weights
    invV         = auxdata{4};
     
    % Extract calibrated parameters from scf
    mhat_scf_param4         = auxdata{5};    

    
    % Number of moments
    n           = length(mhat);
    
    % Number of parameters
    k           = length(theta);

    % Extract individual parameters from parameter vector
    pi_i1       = theta(1);
    pi_i2       = theta(2);
    pi_c1       = theta(3);
    pi_c2       = theta(4);
    
    % Extract individual sample moments from SCF
    mu_r1       = mhat_scf_param4(1);
    mu_r2       = mhat_scf_param4(2);
    mu_a1       = mhat_scf_param4(3);
    mu_a2       = mhat_scf_param4(4);
    var_r1      = mhat_scf_param4(5);
    var_r2      = mhat_scf_param4(6);
    var_a1      = mhat_scf_param4(7);
    var_a2      = mhat_scf_param4(8);
    
    % Extract individual sample moments from SCF
    cov_r1r2    = mhat_scf(1);
    cov_r1a1    = mhat_scf(2);
    cov_r1a2    = mhat_scf(3);
    cov_r2a1    = mhat_scf(4);
    cov_r2a2    = mhat_scf(5);
    cov_a1a2    = mhat_scf(6);

    % Model moments
    mu_y1       = mu_r1+mu_a1;
    mu_y2       = mu_r2+mu_a2;
    mu_a        = sa1*mu_a1+(1-sa1)*mu_a2;
    mu_ri       = pi_i1*mu_r1+pi_i2*mu_r2;
    mu_rc       = pi_c1*mu_r1+pi_c2*mu_r2;
    var_y1      = var_r1+var_a1+2*cov_r1a1;
    var_y2      = var_r2+var_a2+2*cov_r2a2;
    var_a       = (sa1^2)*var_a1+((1-sa1)^2)*var_a2+2*sa1*(1-sa1)*cov_a1a2;
    var_ri      = (pi_i1^2)*var_r1+(pi_i2^2)*var_r2+2*pi_i1*pi_i2*cov_r1r2;
    var_rc      = (pi_c1^2)*var_r1+(pi_c2^2)*var_r2+2*pi_c1*pi_c2*cov_r1r2;
    cov_y1y2    = cov_r1r2+cov_r1a2+cov_r2a1+cov_a1a2;
    cov_y1a     = sa1*cov_r1a1+sa1*var_a1+...
                (1-sa1)*cov_r1a2+(1-sa1)*cov_a1a2;
    cov_y1ri    = pi_i1*var_r1+pi_i1*cov_r1a1+...
                pi_i2*cov_r1r2+pi_i2*cov_r2a1;
    cov_y1rc    = pi_c1*var_r1+pi_c1*cov_r1a1+...
                pi_c2*cov_r1r2+pi_c2*cov_r2a1;
    cov_y2a     = sa1*cov_r2a1+sa1*cov_a1a2+...
                (1-sa1)*cov_r2a2+(1-sa1)*var_a2;
    cov_y2ri    = pi_i1*cov_r1r2+pi_i1*cov_r1a2+...
                pi_i2*var_r2+pi_i2*cov_r2a2;
    cov_y2rc    = pi_c1*cov_r1r2+pi_c1*cov_r1a2+...
                pi_c2*var_r2+pi_c2*cov_r2a2;
    cov_ari     = pi_i1*sa1*cov_r1a1+pi_i1*(1-sa1)*cov_r1a2+...
                pi_i2*sa1*cov_r2a1+pi_i2*(1-sa1)*cov_r2a2;
    cov_arc     = pi_c1*sa1*cov_r1a1+pi_c1*(1-sa1)*cov_r1a2+...
                pi_c2*sa1*cov_r2a1+pi_c2*(1-sa1)*cov_r2a2;
    cov_rirc    = pi_i1*pi_c1*var_r1+(pi_i1*pi_c2+pi_i2*pi_c1)*cov_r1r2+...
                pi_i2*pi_c2*var_r2;

    % Vector of model moments
    m           = [mu_y1;mu_y2;mu_a;mu_ri;mu_rc;var_y1;cov_y1y2;...
                cov_y1a;cov_y1ri;cov_y1rc;var_y2;cov_y2a;cov_y2ri;...
                cov_y2rc;var_a;cov_ari;cov_arc;var_ri;cov_rirc;var_rc];
            
    % Difference between model and sample moments
    d           = m-mhat;
    
    % Matrix of first derivatives of model moments wrt model parameters
    G           = zeros(n,k);
%     G(1,1)      = 1;
%     G(1,3)      = 1;
%     G(2,2)      = 1;
%     G(2,4)      = 1;
%     G(3,3)      = sa1;
%     G(3,4)      = 1-sa1;
%     G(4,1)      = pi_i1;
%     G(4,2)      = pi_i2;
    G(4,1)      = mu_r1;
    G(4,2)     = mu_r2;
%     G(5,1)      = pi_c1;
%     G(5,2)      = pi_c2;
    G(5,3)     = mu_r1;
    G(5,3)     = mu_r2;
%     G(6,5)      = 1;
%     G(6,7)      = 1;
%     G(8,7)      = sa1;
%     G(9,5)      = pi_i1;
    G(9,1)      = var_r1+cov_r1a1;
    G(9,2)     = cov_r1r2+cov_r2a1;
%     G(10,5)     = pi_c1;
    G(10,3)    = var_r1+cov_r1a1;
    G(10,4)    = cov_r1r2+cov_r2a1;
%     G(11,6)     = 1;
%     G(11,8)     = 1;
%     G(12,8)     = 1-sa1;
%     G(13,6)     = pi_i2;
    G(13,1)     = cov_r1r2+cov_r1a2;
    G(13,2)    = var_r2+cov_r2a2;
%     G(14,6)     = pi_c2;
    G(14,3)    = cov_r1r2+cov_r1a2;
    G(14,4)    = var_r2+cov_r2a2;
%     G(15,7)     = sa1^2;
%     G(15,8)     = (1-sa1)^2;
    G(16,1)     = sa1*cov_r1a1+(1-sa1)*cov_r1a2;
    G(16,2)    = sa1*cov_r2a1+(1-sa1)*cov_r2a2;
    G(17,3)    = sa1*cov_r1a1+(1-sa1)*cov_r1a2;
    G(17,4)    = sa1*cov_r2a1+(1-sa1)*cov_r2a2;
%     G(18,5)     = pi_i1^2;
%     G(18,6)     = pi_i2^2;
    G(18,1)     = 2*pi_i1*var_r1+2*pi_i2*cov_r1r2;
    G(18,2)    = 2*pi_i2*var_r2+2*pi_i1*cov_r1r2;
%     G(19,5)     = pi_i1*pi_c1;
%     G(19,6)     = pi_i2*pi_c2;
    G(19,1)     = pi_c1*var_r1+pi_c2*cov_r1r2;
    G(19,2)    = pi_c2*var_r2+pi_c1*cov_r1r2;
    G(19,3)    = pi_i1*var_r1+pi_i2*cov_r1r2;
    G(19,4)    = pi_i2*var_r2+pi_i1*cov_r1r2;
%     G(20,5)     = pi_c1^2;
%     G(20,6)     = pi_c2^2;
    G(20,3)    = 2*pi_c1*var_r1+2*pi_c2*cov_r1r2;
    G(20,4)    = 2*pi_c2*var_r2+2*pi_c1*cov_r1r2;
    
    % Gradient
    g           = 2*(d')*G;
    